Tägliche Wetterdaten

Um mit den täglichen Wetterdaten arbeiten zu können, muss zunächst die Spalte date in das richtige Datumsformat überführt werden. Die Temperaturangaben in Fahrenheit werden in Celsius, die Menge des Regens und des Schnees in Inches in Millimeter umgerechnet. Kann die Menge des gefallenen Regens oder Schnees nicht gemessen werden, wird dies durch ein T gekennzeichnet. Um mit den Daten problemlos arbeiten zu können, wird das T durch den Wert 0.01 ersetzt. So ist ein numerischer Wert vorhanden, der kennzeichnet, dass Regen oder Schnee gefallen ist. Jedoch ist dieser sehr gering, sodass er keinen unverhältnismäßigen Einfluss auf etwaige Analyseergebnisse hat. Die T-Werte werden zunächst durch den Wert 100 ersetzt, da so Probleme beim Umrechnen umgangen werden, die durch z.B. NAs entstehen. In keiner der Spalten kommt vorher der Wert 100 vor, sodass nach dem Umrechnen weiterhin klar ist, welche Werte vorher ein T waren. Nach dem Umrechnen haben diese den Wert 2540 und können problemlos mit 0.01 ersetzt werden.

#Datum richtig formatieren
weather$date <- dmy(weather$date)

#Temperatur in Celsius umrechnen
weather <- weather %>%
  mutate(maximum.temperature = round((maximum.temperature-32)*5/9, 2),
         minimum.temperature = round((minimum.temperature-32)*5/9, 2),
         average.temperature = round((average.temperature-32)*5/9, 2))

levels(weather$precipitation)
levels(weather$snow.fall)
levels(weather$snow.depth)

#Niederschlagsmenge formatieren, T-Werte durch 100 ersetzen und zurückformatieren
weather$precipitation <- as.character(weather$precipitation)
weather$precipitation[weather$precipitation == "T"] <- 100
weather$precipitation <- as.double(weather$precipitation)

#Schneefall und -tiefe formatieren, T-Werte durch 100 ersetzten und zurückformatieren
weather$snow.fall <- as.character(weather$snow.fall)
weather$snow.fall[weather$snow.fall == "T"] <- 100
weather$snow.fall <- as.double(weather$snow.fall)

weather$snow.depth <- as.character(weather$snow.depth)
weather$snow.depth[weather$snow.depth == "T"] <- 100
weather$snow.depth <- as.double(weather$snow.depth)

#Regen und Schnee in mm umrechen
weather <- weather %>%
  mutate(precipitation = round(precipitation*25.4, 2),
         snow.fall = round(snow.fall*25.4, 2),
         snow.depth = round(snow.depth*25.4, 2))

#Die ehemaligen T-Werte (jetzt 2540) mit 0.01 ersetzten
max(weather$precipitation)
weather$precipitation[weather$precipitation == 2540] <- 0.01

max(weather$snow.fall)
weather$snow.fall[weather$snow.fall == 2540] <- 0.01

max(weather$snow.depth)
weather$snow.depth[weather$snow.depth == 2540] <- 0.01

# Wetter gespeichert

saveRDS(weather, "Data/weather_daily_2016.rds")

weather <- readRDS("Data/weather_daily_2016.rds")

Die Wetterdaten für das Jahr 2017, werden genauso aufbereitet wie die vorherigen. Die Maßeinheiten Fahrenheit und Inches werden wieder in Celsius und Millimeter umgerechnet. Da die Spalte mit der täglichen Durchschnittstemperatur ausschließlich fehlende Werte enthält, wird sie aus dem Mittelwert der minimal und maximal Temperatur ermittelt. Dies ist ebenfalls das Vorgehen des NCDC. Die Spalten Niederschlagsmenge, Schneefall und Schneetiefe enthalten im Jahr 2017 keine T-Werte wie noch im Jahr 2016 und müssen hier nicht dementsprechend aufbereitet werden.

#Daten laden
weather17 <- read_csv("Data/daily_weather_2017.csv")

#Daten für 2017 auswählen, Spaltennamen dem Vorjahr anpassen und Durchschnittstemperatur berechnen
weather17 <- weather17 %>%
  filter(year(DATE) == 2017) %>%
  select(date = DATE,
         maximum.temperature = TMAX,
         minimum.temperature = TMIN,
         precipitation = PRCP,
         snow.fall = SNOW,
         snow.depth = SNWD) %>%
  mutate(average.temperature = (minimum.temperature+maximum.temperature)/2)

#Fehlende Werte überprüfen
apply(df, 2, function(x) sum(is.na(x)))

#Einheiten umrechnen
weather17 <- weather17 %>%
  mutate(maximum.temperature = round((maximum.temperature-32)*5/9, 2),
         minimum.temperature = round((minimum.temperature-32)*5/9, 2),
         average.temperature = round((average.temperature-32)*5/9, 2),
         precipitation = round(precipitation*25.4, 2),
         snow.fall = round(snow.fall*25.4, 2),
         snow.depth = round(snow.depth*25.4, 2))

#Data Frame wie Vorjahr sortieren
weather17 <- weather17[,c(1:3,7,4:6)]

#Datensatz speichern
saveRDS(weather17, "Data/weather_daily_2017.RDS")

weather_2017 <- read_rds("Data/weather_daily_2017.RDS")

Stündliche Wetterdaten

Auch die stündlichen Wetterdaten müssen vor Gebrauch bearbeitet werden. Der Datensatz hat in vielen Spalten sehr viele fehlende Werte. In einem ersten Schritt werden die Spalten ausgesucht, die weiterhin betrachtet werden sollen. Die Spalten Datum, Temperatur in Grad Celsius, Windgeschwindigkeit in Km/h, sowie die binären Spalten Nebel, Regen, Schnee, Hagel, Gewitter und Tornado werden weiterhin genutzt. Die Spalte, die z.B den gefallenen Regen in mm enthält, weißt mehr fehlende Werte auf, als Werte da sind. Da dadurch ein sinnvolles interpolieren der Werte nicht mehr möglich ist, wird diese Variable, sowie einige andere von der weiteren Betrachtung ausgeschlossen. Die Spalten Temperatur und Windgeschwindigkeit weisen einige fehlende Werte auf, die aber aufgrund ihrer geringen Anzahl sinnvoll interpoiert werden können und werden. Hierfür wird die Spline-Interpolation genutzt. Auch in diesem Datensatz wird die Spalte, die Datum und Uhrzeit enthält aufgespalten, sodass Datum und Stunde als getrennte Spalten übrig bleiben. Im letzten Schritt der Aufbereitung werden doppelte Beobachtungen behandelt. Für manche Stunden sind verschiedene Temperaturen und Windgeschwindigkeiten bekannt, die sich jeweils nur marginal unterscheiden. Es wird jeweils der erste gemessene Wert für die weiteren Analysen beibehalten.

#Nur Wetterdaten aus 2016 behalten
hourly_weather <- hourly_weather %>%
  filter(date(pickup_datetime) > "2015-12-31" & date(pickup_datetime) < "2017-01-01")

#Fehlende Werte überprüfen
apply(hourly_weather, 2, function(x) sum(is.na(x)))

#Spalten auswählen
hourly_weather <- hourly_weather %>%
  select(pickup_datetime, tempm, wspdm, fog, rain, snow, hail, thunder, tornado)

#Auf fehlende Werte überprüfen und diese Interpolieren
apply(hourly_weather, 2, function(x) sum(is.na(x)))
hourly_weather[is.na(hourly_weather$tempm),]
hourly_weather$tempm <- na_interpolation(hourly_weather$tempm, option = "spline")
hourly_weather[is.na(hourly_weather$wspdm),]
hourly_weather$wspdm <- na_interpolation(hourly_weather$wspdm, option = "spline")
apply(hourly_weather, 2, function(x) sum(is.na(x)))

#Spalten umbenennen
hourly_weather <- hourly_weather %>%
  rename(datetime = pickup_datetime,
         temp = tempm,
         wspd = wspdm)

#Datum und Stunde trennen
hourly_weather <- hourly_weather %>%
  mutate(date = as_date(datetime),
         hour = hour(datetime))
hourly_weather$datetime <- NULL         

#Reihenfolge der Spalten festlegen und nur einen Messwert je Stunde behalten
hourly_weather <- hourly_weather[,c(9,10,1:8)]
hourly_weather <- hourly_weather %>% distinct(date, hour, .keep_all = T)

#Daten speichern
saveRDS(hourly_weather, "hourly_weather.RDS")

Aufgrund der dünnen Datenlage zum stündlichen Wetter im Jahr 2017, werden zwei andere Datensätze für die Jahre 2016 und 2017 genutzt, die das stündliche Wetter enthalten.

Auch diese müssen zunächst aufbereitet werden. Da sehr viele Variablen vorhanden sind, die zum Teil viele fehlende Werte enthalten werden die stündlichen Daten für die Temperatur an der Messstation in Fahrenheit, die Windgeschwindigkeit in Meilen pro Stunde, sowie die Niederschlagsmenge in inches ausgewählt. Fahrenheit wird in Celsius, Meilen pro Stunde in Kilometer pro Stunde und inches in Millimeter umgerechnet. Kann die Menge des gefallenen Regens nicht gemessen werden, wird dies durch ein T gekennzeichnet. Um mit den Daten problemlos arbeiten zu können, wird das T durch den Wert 0.01 ersetzt. So ist ein numerischer Wert vorhanden, der kennzeichnet, dass Regen gefallen ist. Jedoch ist dieser sehr gering, sodass er keinen unverhältnismäßigen Einfluss auf etwaige Analyseergebnisse hat. Die T-Werte werden zunächst durch den Wert 100 ersetzt, da so Probleme beim Umrechnen umgangen werden, die durch z.B. NAs entstehen. In keiner der Spalten kommt vorher der Wert 100 vor, sodass nach dem Umrechnen weiterhin klar ist, welche Werte vorher ein T waren. Nach dem Umrechnen haben diese den Wert 2540 und können problemlos mit 0.01 ersetzt werden. Sind mehrere Werte für eine Stunde vorhanden, so wird nur der erste Wert weiter berücksichtigt. Verbliebene fehlende Werte werden mit Hilfe der Spline-Interpolation ersetzt. Dies geschieht für die Jahre 2016 und 2017 getrennt.

#Daten einlesen
df <- read_csv("Data/hourly_weather_new.csv")

#stündliche Daten auswählen
df <- df %>% select(DATE, HourlyAltimeterSetting, HourlyDewPointTemperature, HourlyDryBulbTemperature, 
                    HourlyPrecipitation, HourlyPresentWeatherType, HourlyPressureChange,   HourlyPressureTendency,
                    HourlyRelativeHumidity, HourlySeaLevelPressure, HourlySkyConditions,  HourlyStationPressure,
                    HourlyVisibility, HourlyWetBulbTemperature, HourlyWindDirection, HourlyWindGustSpeed, 
                    HourlyWindSpeed)

#Fehlende Werte überprüfen
apply(df, 2, function(x) sum(is.na(x)))

#Temperatur, Windgeschwindigkeit und Niederschlagsmenge auswählen
df <- df %>% select(DATE, HourlyDryBulbTemperature, HourlyPrecipitation, HourlyWindSpeed)

#Datensatz nach Jahren 2016 und 2017 aufsplitten
df16 <- df %>%
  filter(year(DATE)==2016)
df17 <- df %>%
  filter(year(DATE)==2017)

#Jeweils fehlende Werte überprüfen
apply(df16, 2, function(x) sum(is.na(x)))
apply(df17, 2, function(x) sum(is.na(x)))

#2016
#Datum und Stunde trennen, sowie Uhrzeit als eigene Spalte behalten
df16 <- df16 %>%
  mutate(date = as_date(DATE),
         hour = hour(DATE),
         time = paste(hour(DATE), minute(DATE), second(DATE), sep = ":"))

#Alle Zeilen mit Uhrzeit 23:59:00 entfernen, da letzter Tagesreport mit allen Werten NA
df16 <- df16 %>% 
  filter(time != "23:59:0")
df16$time <- NULL
df16$DATE <- NULL

#Temperatur und Windgeschwindigkeit umrechnen
df16 <- df16 %>%
  mutate(temp = round((HourlyDryBulbTemperature-32)*5/9, 2),
         wdsp = round(HourlyWindSpeed*1.609344, 2))
df16$HourlyDryBulbTemperature <- NULL
df16$HourlyWindSpeed <- NULL

#Niederschlagsmenge umbenennen, Werte mit Buchstaben korrigieren, T durch 1000 ersetzten, in mm umrechnen und die ehemaligen T-Werte in 0.01mm ändern
df16 <- df16 %>%
  rename(precip = HourlyPrecipitation)

df16$precip <- str_replace_all(df16$precip, c("s" = ""))
df16$precip[df16$precip == "T"] <- 1000
df16$precip <- as.double(df16$precip)

df16 <- df16 %>%
  mutate(precip = round(precip*25.4, 2))
max(df16$precip, na.rm = T)
df16$precip[df16$precip == 25400] <- 0.01

#Reihenfolge der Spalten ändern und nur eine Zeile pro Stunde behalten
df16 <- df16[,c(2,3,4,5,1)]
df16 <- df16 %>% 
  distinct(date, hour, .keep_all = T)

#Fehlende Werte interpolieren
df16[,3:5] <- na_interpolation(df16[,3:5], option = "spline")
apply(df16, 2, function(x) sum(is.na(x)))

#2017
#Datum und Stunde trennen, sowie Uhrzeit als eigene Spalte behalten
df17 <- df17 %>%
  mutate(date = as_date(DATE),
         hour = hour(DATE),
         time = paste(hour(DATE), minute(DATE), second(DATE), sep = ":"))

#Alle Zeilen mit Uhrzeit 23:59:00 entfernen, da letzter Tagesreport mit allen Werten NA
df17 <- df17 %>% 
  filter(time != "23:59:0")
df17$time <- NULL
df17$DATE <- NULL

#Temperatur und Windgeschwindigkeit umrechnen
df17 <- df17 %>%
  mutate(temp = round((HourlyDryBulbTemperature-32)*5/9, 2),
         wdsp = round(HourlyWindSpeed*1.609344, 2))
df17$HourlyDryBulbTemperature <- NULL
df17$HourlyWindSpeed <- NULL

#Niederschlagsmenge umbenennen, Werte mit Buchstaben korrigieren, T durch 1000 ersetzten, in mm umrechnen und die ehemaligen T-Werte in 0.01mm ändern
df17 <- df17 %>%
  rename(precip = HourlyPrecipitation)

df17$precip <- str_replace_all(df17$precip, c("s" = ""))
df17$precip[df17$precip == "T"] <- 1000
df17$precip <- as.double(df17$precip)

df17 <- df17 %>%
  mutate(precip = round(precip*25.4, 2))
max(df17$precip, na.rm = T)
df17$precip[df17$precip == 25400] <- 0.01

#Reihenfolge der Spalten ändern und nur eine Zeile pro Stunde behalten
df17 <- df17[,c(2,3,4,5,1)]
df17 <- df17 %>% 
  distinct(date, hour, .keep_all = T)

#Fehlende Werte interpolieren
df17[,3:5] <- na_interpolation(df17[,3:5], option = "spline")
apply(df17, 2, function(x) sum(is.na(x)))

#Daten speichern
saveRDS(df16, "Data/hourly_weather_2016.RDS")
saveRDS(df17, "Data/hourly_weather_2017.RDS")

hourly_weather_16 <- read_rds("Data/hourly_weather_2016.RDS")
hourly_weather_17 <- read_rds("Data/hourly_weather_2017.RDS")

Im weiteren Verlauf dieser Arbeit wird mit den Datensätzen hourly_weather_16 und hourly_weather_2017 gearbeitet.

Stündliche Daten

Um die erste Forschungsfrage beantworten zu können, müssen die Daten zunächst so aggregiert werden, dass die Anzahl der Fahrten pro Stunde und Station sichtbar werden. Hier werden zwei Datensätze erstellt: Der Erste hält die Anzahl der pro Stunde losgefahrenen Radfahrer je Station fest (hourly_starts). Der Zweite hält die Anzahl der pro Stunde ankommenden Radfahrer je Station fest (hourly_stops). Die Anzahl der Fahrer wird mit der Variable user_count beschrieben. Die Variablen avg_age und avg_tripduration beschreiben jeweils das durchschnittliche Alter in Jahren, sowie die durchschnittliche Fahrtdauer in Sekunden. Die Variable weekend gibt an, ob der jeweilige Tag ein Wochenendtag (Samstag oder Sonntag) ist (1 = ja, 0 = nein). In den Variablen male_user_count, female_user_count und undefined_user_count geben die Anzahl der Fahrer nach Geschlecht gesplittet an. subscriber_count und customer_count geben an, wie viele der Fahrer ein Abo hatten bzw. einen 1-Tages oder 3-Tages-Pass. In den letzten 5 genannten Variablen werden fehlende Werte durch 0 ersetzt, da sich diese durch das Pivotieren der jeweiligen Data Frames ergeben und für die Anzahl der Fahrer stehen.

#Datensatz mit getrenntem Datum und Stunde erstellen
bike_q1 <- bike
bike_q1 <- bike_q1 %>%
  mutate(start_date = as_date(starttime),
         start_hour = hour(starttime),
         stop_date = as_date(stoptime),
         stop_hour = hour(stoptime))

#Datensatz für stündliche Starts je Station erstellen und Spalten für die Anzahl der losfahrenden Radfahrer, deren durchschnittliches Alter und die durchschnittliche Fahrtdauer einfügen, sowie Variable für Wochenende
hourly_starts <- bike_q1 %>%
  group_by(start_date, start_hour, start_station_id, start_station_name, start_station_latitude,
           start_station_longitude) %>%
  summarise(user_count = n(),
            avg_age = round(2016 - mean(birth_year)),
            avg_tripduration = mean(tripduration)) %>%
  mutate(weekend = ifelse(wday(start_date) == 6 | wday(start_date) == 7, 1,0))

#Datensatz erstellen der nach Stunde, Station und Geschlecht getrennt Fahrer zählt
gender <- bike_q1 %>%
  group_by(start_date, start_hour, start_station_id, gender) %>%
  summarise(user = n())
gender <- pivot_wider(gender, names_from = gender, values_from = user)

#Beide Datensätze verbinden und redundante Spalten löschen
hourly_starts <- cbind(hourly_starts, gender)
hourly_starts$start_date1 <- NULL
hourly_starts$start_hour1 <- NULL
hourly_starts$start_station_id1 <- NULL

#Datensatz erstellen, der nach Stunde, Station und Nutzertyp getrennt Fahrer zählt
usertype <- bike_q1 %>%
  group_by(start_date, start_hour, start_station_id, usertype) %>%
  summarise(user = n())
usertype <- pivot_wider(usertype, names_from = usertype, values_from = user)

#Beide Datensätze verbinden und redundante Spalten löschen
hourly_starts <- cbind(hourly_starts, usertype)
hourly_starts$start_date1 <- NULL
hourly_starts$start_hour1 <- NULL
hourly_starts$start_station_id1 <- NULL

#Spalten umbenennen
hourly_starts <- hourly_starts %>%
  rename(male_user_count = "1",
         female_user_count = "2",
         undefined_user_count = "0",
         subscriber_count = Subscriber,
         customer_count = Customer)

#Fehlende Werte ersetzen
hourly_starts$male_user_count[is.na(hourly_starts$male_user_count)] <- 0
hourly_starts$female_user_count[is.na(hourly_starts$female_user_count)] <- 0
hourly_starts$undefined_user_count[is.na(hourly_starts$undefined_user_count)] <- 0
hourly_starts$subscriber_count[is.na(hourly_starts$subscriber_count)] <- 0
hourly_starts$customer_count[is.na(hourly_starts$customer_count)] <- 0

#nicht mehr benötigte Datensätze löschen
rm(gender)
rm(usertype)

#Datensatz für stündliche Stops je Station erstellen und Spalten für die Anzahl der losfahrenden Radfahrer, deren durchschnittliches Alter und die durchschnittliche Fahrtdauer einfügen, sowie Variable für Wochenende
hourly_stops <- bike_q1 %>%
  group_by(stop_date, stop_hour, end_station_id, end_station_name, end_station_latitude,
           end_station_longitude) %>%
  summarise(user_count = n(),
            avg_age = round(2016 - mean(birth_year)),
            avg_tripduration = mean(tripduration)) %>%
  mutate(weekend = ifelse(wday(stop_date) == 6 | wday(stop_date) == 7, 1,0))

#Datensatz erstellen der nach Stunde, Station und Geschlecht getrennt Fahrer zählt
gender <- bike_q1 %>%
  group_by(stop_date, stop_hour, end_station_id, end_station_name, end_station_latitude,
           end_station_longitude,gender) %>%
  summarise(user = n())
gender <- pivot_wider(gender, names_from = gender, values_from = user)

#Beide Datensätze verbinden und redundante Spalten löschen
hourly_stops <- cbind(hourly_stops, gender)
hourly_stops$stop_date1 <- NULL
hourly_stops$stop_hour1 <- NULL
hourly_stops$end_station_id1 <- NULL
hourly_stops$end_station_name1 <- NULL 
hourly_stops$end_station_latitude1 <- NULL
hourly_stops$end_station_longitude1 <- NULL

#Datensatz erstellen, der nach Stunde, Station und Nutzertyp getrennt Fahrer zählt
usertype <- bike_q1 %>%
  group_by(stop_date, stop_hour, end_station_id, end_station_name, end_station_latitude,
           end_station_longitude, usertype) %>%
  summarise(user = n())
usertype <- pivot_wider(usertype, names_from = usertype, values_from = user)

#Beide Datensätze verbinden und redundante Spalten löschen
hourly_stops <- cbind(hourly_stops, usertype)
hourly_stops$stop_date1 <- NULL
hourly_stops$stop_hour1 <- NULL
hourly_stops$end_station_id1 <- NULL
hourly_stops$end_station_name1 <- NULL 
hourly_stops$end_station_latitude1 <- NULL
hourly_stops$end_station_longitude1 <- NULL

#Spalten umbenennen
hourly_stops <- hourly_stops %>%
  rename(male_user_count = "1",
         female_user_count = "2",
         undefined_user_count = "0",
         subscriber_count = Subscriber,
         customer_count = Customer)

#Fehlende Werte ersetzen
hourly_stops$male_user_count[is.na(hourly_stops$male_user_count)] <- 0
hourly_stops$female_user_count[is.na(hourly_stops$female_user_count)] <- 0
hourly_stops$undefined_user_count[is.na(hourly_stops$undefined_user_count)] <- 0
hourly_stops$subscriber_count[is.na(hourly_stops$subscriber_count)] <- 0
hourly_stops$customer_count[is.na(hourly_stops$customer_count)] <- 0

#Nicht mehr benötigte Datensätze löschen
rm(gender)
rm(usertype)

#Aggregierte Datensätze speichern
saveRDS(hourly_starts, "Data/bike_q1_starts.rds")
saveRDS(hourly_stops, "Data/bike_q1_stops.rds")

hourly_starts <- readRDS("Data/bike_q1_starts.rds")
hourly_stops <- readRDS("Data/bike_q1_stops.rds")

Tägliche Daten

Um die Daten für die zweite Frage zu aggregieren, müssen die Abfahrten und Ankünfte pro Station pro Tag zusammengefasst werden. Deswegen werde ich als erstes das Datum extrahieren damit ich damit nach dem Tag aggregieren kann.

# Kopiere Datensatz in Question2 Datensatz und wandel Date/time Spalte in Date spalte um und füge Sie hinzu
bike_q2 <- bike
bike_q2 <- bike_q2 %>%
  mutate(date_start=date(bike_q2$starttime),
         date_stop=date(bike_q2$stoptime))


# Extrahiere das Gender aus dem Datensatz für die Starts
bike_q2_starts_gender <- bike_q2 %>%
  group_by(date_start, start_station_name, start_station_latitude, start_station_longitude, gender) %>%
  summarise( anzahl = n() )

# Extrahiere das Gender aus dem Datensatz für die Stops
bike_q2_stops_gender <- bike_q2 %>%
  group_by(date_stop, end_station_name, end_station_latitude, end_station_longitude, gender) %>%
  summarise( anzahl = n() )

# Tidy Gender Tabelle für starts und stops
gender_q2_starts <- pivot_wider(bike_q2_starts_gender, names_from = gender, values_from = anzahl)

gender_q2_stops <- pivot_wider(bike_q2_stops_gender, names_from = gender, values_from = anzahl)


# Aggregiere Starts und Stops nach Datum und Station

bike_q2_starts <- bike_q2 %>%
  group_by(date_start, start_station_name, start_station_latitude, start_station_longitude) %>%
  summarise(avr_age = round(2016-mean(birth_year)), avg_tripduration = mean(tripduration), anzahl = n() )
bike_q2_stops <- bike_q2 %>%
  group_by(date_stop, end_station_name, end_station_latitude, end_station_longitude) %>%
  summarise(avr_age = round(2016-mean(birth_year)), avg_tripduration = mean(tripduration), anzahl = n() )


# Füge den agregierten Tabellen das Gender hinzu

bike_q2_starts_rdy <- left_join(bike_q2_starts, gender_q2_starts, by = c("date_start", "start_station_name" ))
bike_q2_stops_rdy <- left_join(bike_q2_stops, gender_q2_stops, by = c("date_stop", "end_station_name" ))

# Nicht benötigte Spalten gelöscht und NAs in 0 umgewandelt weil wir wissen das diese 0 sind

bike_q2_starts_rdy$start_station_latitude.y <- NULL
bike_q2_starts_rdy$start_station_longitude.y <- NULL

bike_q2_starts_rdy[is.na(bike_q2_starts_rdy$`1`),"1"] <- 0
bike_q2_starts_rdy[is.na(bike_q2_starts_rdy$`2`),"2"] <- 0
bike_q2_starts_rdy[is.na(bike_q2_starts_rdy$`0`),"0"] <- 0

bike_q2_starts_rdy <- rename(bike_q2_starts_rdy, male = "1")
bike_q2_starts_rdy <- rename(bike_q2_starts_rdy, woman = "2")
bike_q2_starts_rdy <- rename(bike_q2_starts_rdy, undefined = "0")


bike_q2_stops_rdy$end_station_latitude.y <- NULL
bike_q2_stops_rdy$end_station_longitude.y <- NULL

bike_q2_stops_rdy[is.na(bike_q2_stops_rdy$`1`),"1"] <- 0
bike_q2_stops_rdy[is.na(bike_q2_stops_rdy$`2`),"2"] <- 0
bike_q2_stops_rdy[is.na(bike_q2_stops_rdy$`0`),"0"] <- 0

bike_q2_stops_rdy <- rename(bike_q2_stops_rdy, male = "1")
bike_q2_stops_rdy <- rename(bike_q2_stops_rdy, woman = "2")
bike_q2_stops_rdy <- rename(bike_q2_stops_rdy, undefined = "0")

rm(bike_q2, bike_q2_starts, bike_q2_starts_gender, bike_q2_stops, bike_q2_stops_gender, gender_q2_starts, gender_q2_stops)


saveRDS(bike_q2_starts_rdy, "Data/bike_q2_starts.rds")
saveRDS(bike_q2_stops_rdy, "Data/bike_q2_stops.rds")

q2_starts <- read_rds("Data/bike_q2_starts.rds")
q2_stops <- read_rds("Data/bike_q2_stops.rds")

Aufbereitung Für Stündliche Analyse

Für die stündliche Analyse werden die Daten wie folgt aufbereitet. Wie zuvor beim Trainingsdatensatz wird die Spalte starttime in die zwei Spalten Datum und Stunde aufgesplittet. Da eine Betrachtung auf Bezirkseben erfolgt, werden die Bezirke zu dem Datensatz hinzugefügt. Für zwei Stationen ist keine Zuordnung zu einem Bezirk möglich, sodass die fehlenden Werte in explizite fehlende Werte umgewandelt werden. Dies hat zur Folge, dass der Bezirk dieser Stationen “Missing” benannt wird. Da in 2017 neue Stationen auch in Bezirken errichtet wurden, die in 2016 noch über keine Stationen verfügten, werden diese “neuen” Bezirke aus den Daten entfernt. Sie können ohne Grundlage aus dem vorherigen Jahr nicht prognostiziert werden. Wie auch die Trainingsdaten werden die Testdaten nach dem Bezirk, dem Datum und der Stunde gruppiert, um die jeweiligen Nutzerzahlen aufsummieren zu können. Die Variablen, die die Ankunftsstationen betreffen, sowie die Nutzerdaten, wie z.B. Alter, Geschlecht und Fahrtdauer, werden nicht weiter betrachtet. Sie werden nicht weiter in die Analyse einfließen, da es das Ziel ist zukünftige Fahrten zu prognostizieren. Es ist dann noch nicht bekannt welches Alter die Fahrer haben werden, wie lang sie fahren oder welches Geschlecht sie haben. Es werden ebenfalls Variablen für den Wochentag und die Feiertage eingefügt, sowie Wetterdaten angefügt.

#Einen Testdatensatz erstellen
hourly_test <- test_data

#Datum und Stunde trennen
hourly_test <- hourly_test %>%
  mutate(start_date = as_date(starttime),
         start_hour = hour(starttime))

#Bezirke für Stationen hinzufügen
hourly_test <- hourly_test %>%
  left_join(stations_2017, by = c("start_station_name"="start_station_name",
                                  "start_station_latitude"="start_station_latitude",
                                  "start_station_longitude"="start_station_longitude"))

#Fehlende Werte überprüfen
apply(hourly_test, 2, function(x) sum(is.na(x)))

#Überprüfen wo der Bezirk fehlt
hourly_test[is.na(hourly_test$neighborhood),]

#Der Bezirk von zwei Stationen ist nicht bekannt, fehlende Werte explizit machen
hourly_test$neighborhood <- fct_explicit_na(hourly_test$neighborhood)

#Nach Bezirk gruppieren und Nutzerzahlen aufsummieren
hourly_test <- hourly_test %>%
  group_by(start_date, start_hour, neighborhood) %>%
  summarise(user_count = n()) %>%
  mutate(weekday = wday(start_date, abbr = F, label = T, week_start = 1),
         holiday = ifelse(isBizday(as.timeDate(start_date), holidayNYSE(2017))==FALSE,1,0),
         month = month(start_date, abbr = F, label = T),
         week = week(start_date))

#Variablen für Wochentag und Monat als ungeordneten Faktor speichern
hourly_test$weekday <- factor(hourly_test$weekday, ordered = FALSE)
hourly_test$month <- factor(hourly_test$month, ordered = FALSE)

#Wetterdaten hinzufügen
hourly_test <- hourly_test %>%
  left_join(hourly_weather_17, by = c("start_date" = "date", "start_hour" = "hour"))

#Fehlende Werte überprüfen
apply(hourly_starts_nh, 2, function(x) sum(is.na(x)))

#Spalte Neighborhood in character umwandeln, um nicht genutzte level zu eliminieren
hourly_test$neighborhood <- as.character(hourly_test$neighborhood)

#Bezirke ausfindig machen, die in 2017 vorkommen, aber nicht in 2016 und diese eliminieren
setdiff(hourly_test$neighborhood, hourly_starts_nh$neighborhood)
hourly_test <- hourly_test %>%
  filter(neighborhood != "Astoria" & neighborhood != "Crown Heights" & neighborhood != "Ditmars Steinway" &
           neighborhood != "Harlem" & neighborhood != "Morningside Heights" & 
           neighborhood != "Prospect Heights" & neighborhood != "Prospect-Lefferts Gardens")

#Als .RDS speichern
saveRDS(hourly_test, "Data/q1_starts_test.RDS")

hourly_test <- read_rds("Data/q1_starts_test.RDS")

Aufbereitung Für Tägliche Analyse

# Dataframe für aggregierte TEsdaten erstellen

test_data_q2 <- test_data

# Nicht benötigte Columns entfernen

test_data_q2$gender <- NULL
test_data_q2$usertype <- NULL


#Datum trennen
test_data_q2 <- test_data_q2 %>%
  mutate(start_date = as_date(starttime))

# stations 2017 umwandeln

stations_2017_test <- stations_2017
stations_2017_test$start_station_name <- as.character(stations_2017_test$start_station_name)
stations_2017_test$neighborhood <- as.character(stations_2017_test$neighborhood)

#Bezirke für Stationen hinzufügen

test_data_q2_test <- test_data_q2 %>%
  left_join(stations_2017, by = c("start_station_name"="start_station_name",
                                  "start_station_latitude"="start_station_latitude",
                                  "start_station_longitude"="start_station_longitude"))


test_data_q2 <- test_data_q2 %>%
  left_join(stations_2017_test, by = c("start_station_name"="start_station_name",
                                  "start_station_latitude"="start_station_latitude",
                                  "start_station_longitude"="start_station_longitude"))



# Points1 benutzen damit man nur die Stations kriegt die auch in den Trainingsdaten vorhanden sind

'points1_test <- points1
points1_test$start_station_name <- as.character(points1_test$start_station_name)
points1_test$neighborhood <- as.character(points1_test$neighborhood)


points1_test <- points1_test %>%
  rename(start_station_latitude = start_station_latitude.x,
         start_station_longitude = start_station_longitude.x)


#Bezirke für Stationen hinzufügen
test_data_q2 <- test_data_q2 %>%
  left_join(points1_test, by = c("start_station_name"="start_station_name",
                                  "start_station_latitude"="start_station_latitude",
                                  "start_station_longitude"="start_station_longitude"))'

#Fehlende Werte überprüfen
#apply(test_data_q2, 2, function(x) sum(is.na(x)))

#Überprüfen wo der Bezirk fehlt
test_data_q2[is.na(test_data_q2$neighborhood),]

#Der Bezirk von zwei Stationen ist nicht bekannt, fehlende Werte explizit machen
test_data_q2$neighborhood <- fct_explicit_na(test_data_q2$neighborhood)




#Spalte Neighborhood in character umwandeln, um nicht genutzte level zu eliminieren
test_data_q2$neighborhood <- as.character(test_data_q2$neighborhood)


test_data_q2$birth_year <- as.numeric(test_data_q2$birth_year)

# Löschen der fehlenden Werte mit Birthyear

test_data_q2 <- na.omit(test_data_q2)

#Nach Bezirk gruppieren und Nutzerzahlen aufsummieren
test_data_q2 <- test_data_q2 %>%
  group_by(start_date, neighborhood) %>%
  summarise(anzahl = n(), avg_age = round(2016-mean(birth_year)), avg_tripduration = mean(tripduration)) %>%
  mutate(weekday = wday(start_date, abbr = F, label = T, week_start = 1),
         holiday = ifelse(isBizday(as.timeDate(start_date), holidayNYSE(2017))==FALSE,1,0),)

#Fehlende Werte überprüfen
apply(test_data_q2, 2, function(x) sum(is.na(x)))

#Wetter hinzufügen

test_data_q2 <- left_join(test_data_q2, weather_2017, by = c("start_date" = "date"))

# Columns wie bei den Trainingsdaten anordnen

test_data_q2 <- test_data_q2[, c(1,2,3,5,4,8,9,10,11,12,13,6,7)]

rm(stations_2017_test, test_data_q2_test, points, points1_test)


#Bezirke ausfindig machen, die in 2017 vorkommen, aber nicht in 2016 und diese eliminieren
setdiff(test_data_q2$neighborhood, q2_starts_w2_agg$neighborhood)
test_data_q2 <- test_data_q2 %>%
  filter(neighborhood != "Astoria" & neighborhood != "Crown Heights" & neighborhood != "Ditmars Steinway" &
           neighborhood != "Harlem" & neighborhood != "Morningside Heights" & 
           neighborhood != "Prospect Heights" & neighborhood != "Prospect-Lefferts Gardens")


#Monat und Woche als Variablen eingefügt


test_data_q2 <- test_data_q2 %>%
  mutate(week = week(start_date))

test_data_q2 <- test_data_q2 %>%
  mutate(month = month(start_date, abbr = F, label = T))


#Variable für Wochentag als ungeordneten Faktor speichern
test_data_q2$weekday <- factor(test_data_q2$weekday, ordered = FALSE)
test_data_q2$month <- factor(test_data_q2$month, ordered = FALSE)

#Als .RDS speichern
saveRDS(test_data_q2, "Data/test_data_q2.RDS")

test_data_q2 <- read_rds("Data/test_data_q2.RDS")

Stündliche Regression

Für die stündlichen Nutzerdaten werden zwei Regressionen durchgeführt. Dies ist der Fall, da sich die Granularität der beiden Datensätze voneinander unterscheidet. Für einen Datensatz sollen die Nutzerzahlen je Bezirk bestimmt werden, für den anderen Datensatz die Nutzerzahlen je Station.

In einem ersten Schritt wird eine lineare Regression auf Bezirksebene durchgeführt. Die Variable start_date wird hierbei nicht berücksichtigt, da sich mit der Variable week eine Kollinearität ergibt. Diese wird weiterhin berücksichtigt.

#Regressionsmodell mit Bezirken
model_q1_nh <- as.formula(user_count ~ . -start_date)
lin_reg_nh <- lm(model_q1_nh, hourly_starts_nh)
summary(lin_reg_nh)

## 
## Call:
## lm(formula = model_q1_nh, data = hourly_starts_nh)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -178.59  -20.50   -3.77   13.40 1211.36 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    -32.377044  29.950792  -1.081 0.279694    
## start_hour                       1.596521   0.014151 112.824  < 2e-16 ***
## neighborhoodBattery Park City   30.574609  29.951256   1.021 0.307344    
## neighborhoodBedford-Stuyvesant  14.798440  29.951009   0.494 0.621244    
## neighborhoodBoerum Hill          1.074560  29.952204   0.036 0.971381    
## neighborhoodBrooklyn Heights    10.532346  29.951495   0.352 0.725103    
## neighborhoodCarroll Gardens      3.680145  29.961980   0.123 0.902244    
## neighborhoodCentral Park        21.123775  29.951437   0.705 0.480644    
## neighborhoodChelsea            161.927113  29.950651   5.406 6.43e-08 ***
## neighborhoodChinatown           20.242432  29.950969   0.676 0.499135    
## neighborhoodCivic Center         9.095917  29.951579   0.304 0.761366    
## neighborhoodClinton Hill        10.964039  29.951166   0.366 0.714318    
## neighborhoodCobble Hill         -4.640456  29.964079  -0.155 0.876926    
## neighborhoodColumbia St         -5.104182  29.955208  -0.170 0.864701    
## neighborhoodDowntown Brooklyn    7.711770  29.951435   0.257 0.796812    
## neighborhoodDUMBO               -0.755089  29.952864  -0.025 0.979888    
## neighborhoodEast Harlem          2.876974  29.959715   0.096 0.923498    
## neighborhoodEast Village        98.643901  29.950655   3.294 0.000989 ***
## neighborhoodFinancial District  43.164296  29.950903   1.441 0.149538    
## neighborhoodFlatiron District   29.692186  29.951166   0.991 0.321514    
## neighborhoodFort Greene         19.213501  29.951149   0.641 0.521202    
## neighborhoodGovernors Island    -9.084588  30.005810  -0.303 0.762072    
## neighborhoodGowanus             -0.832403  29.962571  -0.028 0.977836    
## neighborhoodGramercy            35.363007  29.951008   1.181 0.237725    
## neighborhoodGreenpoint          20.470390  29.951050   0.683 0.494316    
## neighborhoodGreenwich Village   54.927881  29.950799   1.834 0.066664 .  
## neighborhoodHell's Kitchen      79.421283  29.950710   2.652 0.008008 ** 
## neighborhoodKips Bay            40.798231  29.950872   1.362 0.173145    
## neighborhoodLong Island City     7.735734  29.951459   0.258 0.796194    
## neighborhoodLower East Side     53.245860  29.950665   1.778 0.075440 .  
## neighborhoodMidtown            161.697625  29.950678   5.399 6.71e-08 ***
## neighborhoodMurray Hill         38.503191  29.951046   1.286 0.198605    
## neighborhoodNavy Yard           -8.169031  29.956632  -0.273 0.785088    
## neighborhoodNoHo                 4.216448  29.951992   0.141 0.888049    
## neighborhoodNolita               7.350898  29.951643   0.245 0.806127    
## neighborhoodPark Slope           4.490135  29.952755   0.150 0.880838    
## neighborhoodProspect Park       -1.590645  29.965041  -0.053 0.957666    
## neighborhoodRed Hook            -5.768641  29.966483  -0.193 0.847348    
## neighborhoodSoHo                45.179488  29.950870   1.508 0.131440    
## neighborhoodSouth Slope         -5.713108  29.965118  -0.191 0.848793    
## neighborhoodStuyvesant Town     25.498741  29.950964   0.851 0.394576    
## neighborhoodSunset Park        -17.889665  30.026918  -0.596 0.551318    
## neighborhoodTheater District    25.925479  29.951063   0.866 0.386713    
## neighborhoodTribeca             40.076475  29.950927   1.338 0.180874    
## neighborhoodTwo Bridges          1.320849  29.952261   0.044 0.964826    
## neighborhoodUpper East Side     65.114306  29.950815   2.174 0.029703 *  
## neighborhoodUpper West Side     67.902385  29.950825   2.267 0.023383 *  
## neighborhoodVinegar Hill        -0.843232  29.952682  -0.028 0.977541    
## neighborhoodWest Village        64.215558  29.950737   2.144 0.032030 *  
## neighborhoodWilliamsburg        47.658821  29.950780   1.591 0.111557    
## weekdayDienstag                  1.694371   0.356513   4.753 2.01e-06 ***
## weekdayMittwoch                  2.391670   0.355035   6.736 1.63e-11 ***
## weekdayDonnerstag                1.198441   0.353616   3.389 0.000701 ***
## weekdayFreitag                  -1.246903   0.350134  -3.561 0.000369 ***
## weekdaySamstag                   4.543729   0.666768   6.815 9.47e-12 ***
## weekdaySonntag                   2.687937   0.667564   4.026 5.66e-05 ***
## holiday                        -17.852841   0.639749 -27.906  < 2e-16 ***
## monthFebruar                     1.162120   0.603868   1.924 0.054298 .  
## monthMärz                        4.676641   0.818224   5.716 1.09e-08 ***
## monthApril                       8.238749   1.094478   7.528 5.18e-14 ***
## monthMai                         9.243367   1.402555   6.590 4.39e-11 ***
## monthJuni                       12.874937   1.704095   7.555 4.19e-14 ***
## monthJuli                        7.819132   2.015035   3.880 0.000104 ***
## monthAugust                     14.597768   2.328460   6.269 3.63e-10 ***
## monthSeptember                  29.861510   2.624208  11.379  < 2e-16 ***
## monthOktober                    43.355630   2.931681  14.789  < 2e-16 ***
## monthNovember                   46.226275   3.246522  14.239  < 2e-16 ***
## monthDezember                   49.376731   3.568117  13.838  < 2e-16 ***
## week                            -0.907800   0.073679 -12.321  < 2e-16 ***
## temp                             1.677223   0.019238  87.181  < 2e-16 ***
## wdsp                             0.039288   0.006483   6.060 1.36e-09 ***
## precip                          -7.322706   0.224793 -32.575  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 51.86 on 316527 degrees of freedom
## Multiple R-squared:  0.3998, Adjusted R-squared:  0.3996 
## F-statistic:  2969 on 71 and 316527 DF,  p-value: < 2.2e-16

anova(lin_reg_nh)

## Analysis of Variance Table
## 
## Response: user_count
##                  Df    Sum Sq  Mean Sq  F value    Pr(>F)    
## start_hour        1  23765489 23765489 8835.236 < 2.2e-16 ***
## neighborhood     48 466442661  9717555 3612.671 < 2.2e-16 ***
## weekday           6  15949222  2658204  988.234 < 2.2e-16 ***
## holiday           1   3656889  3656889 1359.512 < 2.2e-16 ***
## month            11  33814466  3074042 1142.829 < 2.2e-16 ***
## week              1    111382   111382   41.408 1.237e-10 ***
## temp              1  20385993 20385993 7578.849 < 2.2e-16 ***
## wdsp              1     83747    83747   31.135 2.409e-08 ***
## precip            1   2854353  2854353 1061.156 < 2.2e-16 ***
## Residuals    316527 851411269     2690                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

train_mse <- mean((lin_reg_nh$residuals)^2)
train_rmse <- sqrt(mean((lin_reg_nh$residuals)^2))
train_mse

## [1] 2689.242

train_rmse

## [1] 51.8579

Die Anova deutet darauf hin, dass der Einfluss aller Variablen auf die Anzahl an Radfahrern statistisch signifikant ist, da der p-Wert äußerst gering ist. Die Funktion lm() hat zur Durchführung der Regression Dummyvariablen für die Variablen neighborhood, den Bezirk, und weekday, den Wochentag, erstellt. Sieht man sich die genaue Zusammenfassung des Modells an, wird sichtbar, dass der Einfluss zweier Bezirke statistisch nicht signifikant ist, d.h. zufällig ist. Sie werden dennoch nicht aus dem Modell entfernt, da sich sonst ein unvollständiges Bild der Situation ergeben würde, denn in diesen Bezirken gibt es auch Leihstationen, die genutzt werden. Sieht man jedoch auf den Wert für R-Quadrat, der sich auf 0.3946 beläuft, wird klar, dass das Modell die Situation nicht in ihrer Vollständigkeit erfassen kann. Dieser Wert ist sehr niedrig und sagt aus, dass nur 39.46% der Varianz in den Daten durch dieses Regressionsmodell erklärt werden können. Wirft man einen Blick auf die Visualisierungen im Kapitel Basic Visualizations wird klar weshalb. Die lineare Regression kann nur in eine Richtung vorhersagen treffen. Z.B. Wenn sich x um 1 erhöht/verringert, erhöht/verringert sich y um z. Es ist jedoch zu sehen, dass die Zusammenhänge in den Daten nicht linearer Natur sind, denn wenn die Temperatur eine bestimmte Marke überschreitet fahren weniger Menschen mit dem Rad. Das Selbe lässt sich bei der Uhrzeit beobachten. Morgens fahren viele Menschen Rad, über Tag weniger und am Nachmittag wieder mehr, bis es in der Nacht wieder weniger wird. Es ist hier also wenig verwunderlich, dass ein lineares Modell den Einfluss der einzelnen Variablen nicht exakt wiedergibt.

#Vorhersage mit Testdaten
lin_reg_nh_predictions <- predict(lin_reg_nh, hourly_test)

test_mse <- mean((hourly_test$user_count - lin_reg_nh_predictions) ^ 2)
test_rmse <- sqrt(mean((hourly_test$user_count - lin_reg_nh_predictions) ^ 2))
test_mse

## [1] 3696.784

test_rmse

## [1] 60.80118

#Bestimmen welche Variablen signifikant sind
test_sig <- summary(lin_reg_nh)$coefficients[-1,4] < 0.05

#Nur signifikante Behalten
sig_x <- names(test_sig)[test_sig == TRUE] 
View(sig_x)
sig_x <- str_replace_all(sig_x, c(" " = "_"))
sig_x <- str_replace_all(sig_x, c("-" = "_"))
sig_x <- str_replace_all(sig_x, c("'" = ""))

#Model mit nur signifikanten Variablen erstellen
sig.model <- as.formula(paste("user_count", paste(sig_x, collapse = " + "), sep = " ~ "))

Nachfolgend wir eine lineare Regression auf Stationsebene durchgeführt.

#Regression mit Stations
model_q1 <- as.formula(user_count ~ .)
lin_reg_q1 <- lm(model_q1, hourly_starts)

Da der Datensatz in bearbeiteter Form mehrere Gigabyte groß ist, kann keine Analyse durchgeführt werden. Nachfolgend wird ausschließlich der Datebsatz auf Bezirksebene betrachtet.

Tägliche Regression

q2_starts_w2_agg <- q2_starts_w2_agg %>%
  rename(
    start_date = date
  )

#Regression mit q2_starts_w
options(max.print=450)
model_q2 <- as.formula(anzahl ~ .-date)
lin_reg_q2 <- lm(model_q2, q2_starts_w)
summary(lin_reg_q2)

#Regression mit q2_starts_w2_agg
options(max.print=450)
model_q2_w2_agg <- as.formula(anzahl ~ .-start_date)
lin_reg_q2_w2_agg <- lm(model_q2_w2_agg, q2_starts_w2_agg)
summary(lin_reg_q2_w2_agg)

## 
## Call:
## lm(formula = model_q2_w2_agg, data = q2_starts_w2_agg)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.2750  -0.8151   0.1902   0.7562  13.6364 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    -6.300e+00  1.393e+00  -4.522 6.18e-06 ***
## neighborhoodBattery Park City   4.495e+00  1.316e+00   3.415 0.000639 ***
## neighborhoodBedford-Stuyvesant  3.393e+01  1.315e+00  25.800  < 2e-16 ***
## neighborhoodBoerum Hill         4.082e+00  1.316e+00   3.103 0.001920 ** 
## neighborhoodBrooklyn Heights    9.304e+00  1.317e+00   7.067 1.65e-12 ***
## neighborhoodCarroll Gardens     1.017e+01  1.325e+00   7.677 1.72e-14 ***
## neighborhoodCentral Park        8.539e+00  1.317e+00   6.486 9.11e-11 ***
## neighborhoodChelsea             2.489e+01  1.316e+00  18.915  < 2e-16 ***
## neighborhoodChinatown           5.692e+00  1.316e+00   4.327 1.52e-05 ***
## neighborhoodCivic Center        3.580e+00  1.316e+00   2.720 0.006531 ** 
## neighborhoodClinton Hill        1.100e+01  1.315e+00   8.361  < 2e-16 ***
## neighborhoodCobble Hill         1.314e+00  1.325e+00   0.992 0.321308    
## neighborhoodColumbia St         1.771e+00  1.316e+00   1.346 0.178251    
## neighborhoodDowntown Brooklyn   8.240e+00  1.315e+00   6.264 3.85e-10 ***
## neighborhoodDUMBO               2.661e+00  1.316e+00   2.023 0.043136 *  
## neighborhoodEast Harlem         1.304e+01  1.323e+00   9.859  < 2e-16 ***
## neighborhoodEast Village        2.190e+01  1.315e+00  16.655  < 2e-16 ***
## neighborhoodFinancial District  1.761e+01  1.316e+00  13.385  < 2e-16 ***
## neighborhoodFlatiron District   3.719e+00  1.316e+00   2.827 0.004706 ** 
## neighborhoodFort Greene         1.782e+01  1.315e+00  13.547  < 2e-16 ***
## neighborhoodGovernors Island    1.357e+00  1.328e+00   1.022 0.306708    
## neighborhoodGowanus             1.018e+01  1.324e+00   7.684 1.63e-14 ***
## neighborhoodGramercy            5.745e+00  1.315e+00   4.368 1.26e-05 ***
## neighborhoodGreenpoint          1.885e+01  1.315e+00  14.335  < 2e-16 ***
## neighborhoodGreenwich Village   1.244e+01  1.316e+00   9.459  < 2e-16 ***
## neighborhoodHell's Kitchen      1.841e+01  1.316e+00  13.990  < 2e-16 ***
## neighborhoodKips Bay            7.892e+00  1.315e+00   6.001 2.00e-09 ***
## neighborhoodLong Island City    1.258e+01  1.315e+00   9.561  < 2e-16 ***
## neighborhoodLower East Side     1.838e+01  1.316e+00  13.971  < 2e-16 ***
## neighborhoodMidtown             3.982e+01  1.316e+00  30.260  < 2e-16 ***
## neighborhoodMurray Hill         5.272e+00  1.316e+00   4.006 6.20e-05 ***
## neighborhoodNavy Yard           2.305e+00  1.316e+00   1.752 0.079767 .  
## neighborhoodNoHo                8.298e-01  1.315e+00   0.631 0.528160    
## neighborhoodNolita              1.617e+00  1.315e+00   1.230 0.218801    
## neighborhoodPark Slope          6.550e+00  1.315e+00   4.980 6.44e-07 ***
## neighborhoodProspect Park       3.061e+00  1.327e+00   2.307 0.021041 *  
## neighborhoodRed Hook            1.018e+01  1.325e+00   7.681 1.67e-14 ***
## neighborhoodSoHo                9.642e+00  1.315e+00   7.330 2.42e-13 ***
## neighborhoodSouth Slope         3.359e+00  1.326e+00   2.533 0.011331 *  
## neighborhoodStuyvesant Town     4.813e+00  1.315e+00   3.659 0.000254 ***
## neighborhoodSunset Park         1.625e+00  1.325e+00   1.226 0.220071    
## neighborhoodTheater District    4.833e+00  1.316e+00   3.672 0.000242 ***
## neighborhoodTribeca             1.245e+01  1.316e+00   9.457  < 2e-16 ***
## neighborhoodTwo Bridges         2.593e+00  1.316e+00   1.971 0.048789 *  
## neighborhoodUpper East Side     2.906e+01  1.316e+00  22.085  < 2e-16 ***
## neighborhoodUpper West Side     2.708e+01  1.317e+00  20.563  < 2e-16 ***
## neighborhoodVinegar Hill        2.991e+00  1.315e+00   2.274 0.022960 *  
## neighborhoodWest Village        1.197e+01  1.316e+00   9.098  < 2e-16 ***
## neighborhoodWilliamsburg        4.251e+01  1.315e+00  32.323  < 2e-16 ***
## avg_tripduration               -1.164e-06  6.641e-06  -0.175 0.860811    
## avg_age                         1.503e-01  1.240e-02  12.121  < 2e-16 ***
## maximum.temperature             2.346e-01  2.632e+00   0.089 0.928972    
## minimum.temperature             2.301e-01  2.632e+00   0.087 0.930322    
## average.temperature            -4.564e-01  5.263e+00  -0.087 0.930903    
## precipitation                  -3.343e-04  2.612e-03  -0.128 0.898172    
## snow.fall                       4.226e-03  3.585e-03   1.179 0.238573    
## snow.depth                     -1.536e-02  1.121e-03 -13.700  < 2e-16 ***
## weekdayDienstag                 2.206e-03  7.106e-02   0.031 0.975237    
## weekdayMittwoch                -6.093e-02  7.081e-02  -0.860 0.389568    
## weekdayDonnerstag              -4.959e-02  7.066e-02  -0.702 0.482846    
## weekdayFreitag                 -5.208e-03  7.012e-02  -0.074 0.940803    
## weekdaySamstag                  2.365e-01  1.309e-01   1.808 0.070683 .  
## weekdaySonntag                  2.235e-01  1.309e-01   1.707 0.087834 .  
## holiday                        -1.679e-01  1.250e-01  -1.343 0.179230    
## week                            2.043e-02  1.501e-02   1.362 0.173322    
## monthFebruar                   -1.329e-01  1.257e-01  -1.058 0.290122    
## monthMärz                      -2.071e-01  1.724e-01  -1.201 0.229753    
## monthApril                     -2.552e-01  2.292e-01  -1.113 0.265565    
## monthMai                       -4.006e-01  2.922e-01  -1.371 0.170368    
## monthJuni                      -4.494e-01  3.569e-01  -1.259 0.207972    
## monthJuli                      -3.845e-01  4.215e-01  -0.912 0.361716    
## monthAugust                     2.523e-02  4.867e-01   0.052 0.958658    
## monthSeptember                  6.812e-01  5.446e-01   1.251 0.211049    
## monthOktober                    1.006e+00  6.048e-01   1.663 0.096307 .  
## monthNovember                   8.040e-01  6.676e-01   1.204 0.228480    
## monthDezember                   7.680e-01  7.337e-01   1.047 0.295228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.265 on 15240 degrees of freedom
## Multiple R-squared:  0.9546, Adjusted R-squared:  0.9544 
## F-statistic:  4275 on 75 and 15240 DF,  p-value: < 2.2e-16

Prediction

#Vorhersage mit Testdaten
lin_reg_q2_predictions <- predict(lin_reg_q2_w2_agg, test_data_q2)

test_q2_mse <- mean((test_data_q2$anzahl - lin_reg_q2_predictions) ^ 2)
test_q2_rmse <- sqrt(mean((test_data_q2$anzahl - lin_reg_q2_predictions) ^ 2))
test_q2_mse

## [1] 1893696

test_q2_rmse

## [1] 1376.116

Stündliche Daten Aufbereiten

Um ein künstliches neuornales Netz zu erhalten, das fehlerfrei läuft, müssen der Trainings- und Testdatensatz der stündlichen Fahrten noch weiter aufbereitet werden. Für die Spalten neighborhood und weekday werden Dummyvariablen erstellt. Dies ist notwendig, da es sich um kategoriale Variablen handelt. Da eine reine Quantifizerung der Variablen durch hohe Zahlenwerte für einige Kategorein eine Verzerrung der Analysen zur Folge hätte, werden binäre Variblen für jede Kategorie erstellt, die die Werte 0 und 1 annehmen. Um bei der späteren Anwendung der künstlichen neuronalen Netze keine Fehler zu produzieren, werden die Dummys effektkodiert. Statt des Wertes 0 wird der Wert -1 genutzt.

#Dummyvariablen erstellen
hourly_starts_nh <- dummy_cols(hourly_starts_nh, c("start_hour", "neighborhood", "weekday", "month", "week"), remove_first_dummy = T, remove_selected_columns = T)
hourly_test <- dummy_cols(hourly_test, c("start_hour", "neighborhood", "weekday", "month", "week"), remove_first_dummy = T, remove_selected_columns = T)

#Namen der Dummys bearbeiten
names(hourly_starts_nh) <- str_replace_all(names(hourly_starts_nh), c(" " = "_"))
names(hourly_test) <- str_replace_all(names(hourly_test), c(" " = "_"))

#Dummys effektkodieren

hourly_starts_nh[,c(3, 7:146)] <- effectcoding(hourly_starts_nh[,c(3, 7:146)])
hourly_test[,c(3, 7:105)] <- effectcoding(hourly_test[,c(4, 7:105)])

Tägliche Daten Aufbereiten

Vorbereiten der des anderen DAtensatzes

# Data Table variante

#library(data.table)

#q2_starts_w_nn <- as.data.table(q2_starts_w)
# Dataframe für LSTM erstellen
q2_starts_w_nn <- q2_starts_w

#Dummyvariablen erstellen
q2_starts_w_nn <- dummy_cols(q2_starts_w, c("weekday", "start_station_name"), remove_first_dummy = T, remove_selected_columns = T)

#Namen der Dummys bearbeiten
names(q2_starts_w_nn) <- str_replace_all(names(q2_starts_w_nn), c(" " = "_"))

#Dummys effektkodieren
#q2_starts_w_nn[,c(11:670)] <- effectcoding(q2_starts_w_nn[,c(11:670)])

Funktioniert nicht!

Neue Variante mit neighborhoods.

Trainingdaten für LSTM vorbereiten

# Data Table variante
#library(data.table)

#q2_starts_w2_nn <- as.data.table(q2_starts_w2_agg)

# Dataframe für LSTM erstellen

q2_starts_w2_nn <- q2_starts_w2_agg

#Dummyvariablen erstellen

q2_starts_w2_nn <- dummy_cols(q2_starts_w2_nn, c("weekday", "neighborhood", "month", "week"), remove_first_dummy = T, remove_selected_columns = T)

#Namen der Dummys bearbeiten
names(q2_starts_w2_nn) <- str_replace_all(names(q2_starts_w2_nn), c(" " = "_"))

#Dummys effektkodieren
q2_starts_w2_nn[,c(11:128)] <- effectcoding(q2_starts_w2_nn[,c(11:128)])

Testdaten für LSTM vorbereiten

# Data Table variante
#test_data_q2_nn <- as.data.table(test_data_q2)

# Dataframe für LSTM erstellen

test_data_q2_nn <- test_data_q2

#Dummyvariablen erstellen
test_data_q2_nn <- dummy_cols(test_data_q2_nn, c("weekday", "neighborhood", "month", "week"), remove_first_dummy = T, remove_selected_columns = T)

#Namen der Dummys bearbeiten
names(test_data_q2_nn) <- str_replace_all(names(test_data_q2_nn), c(" " = "_"))

#Dummys effektkodieren
test_data_q2_nn[,c(11:87)] <- effectcoding(test_data_q2_nn[,c(11:87)])

Vergleich der Datensätze

df1 <- colnames(test_data_q2_nn)
df2 <- colnames(q2_starts_w2_nn)


setdiff(df1,df2)

## character(0)

test_data_q2_nn$neighborhood_Astoria <- NULL
test_data_q2_nn$neighborhood_Crown_Heights <- NULL
test_data_q2_nn$neighborhood_Ditmars_Steinway <- NULL
test_data_q2_nn$neighborhood_Harlem <- NULL
test_data_q2_nn$neighborhood_Morningside_Heights <- NULL
test_data_q2_nn$neighborhood_Prospect_Heights <- NULL
test_data_q2_nn$'neighborhood_Prospect-Lefferts_Gardens' <- NULL


rm(df1,df2)

cross validation